function heat_capacity()
    n = 40;
    p = 0.5;
    b = 0;
    isingmat = isingmat_ini(n, p);
    flip_num = 1;
    T = 5;
    
    %需要调节的参数
    list_len = 15;
    collect_step = 6;
    num = 21;
    gap = 1000;
    img_num = 5;
    
    steps = linspace(0, gap*(num-1), num);
    
    Clist = zeros(1,num);
    
    %en_list = zeros(1,list_len);
    
%     for a = 1:100000
%         [trans_mat, k, m]= trans_state(isingmat, n, T, b, flip_num);
%         isingmat = trans_mat;
%     end 
    mat_list = cell(1,num);
    
    real_gap = gap - list_len * collect_step;

    for i = 1:num
        mat_list{i} = isingmat;
        
        en_list = zeros(1,list_len);
        
        for p = 1:list_len
            for q = 1:collect_step
                [trans_mat, k, m]= trans_state(isingmat, n, T, b, flip_num);
                isingmat = trans_mat;
            end
            en = ising_energy_in(isingmat, n, b);
            en_list(p) = en;
        end
        
        C = var(en_list)/(T*T); 
        Clist(i) = C;
        
        if(i ~= num)
            for j = 1:real_gap
                [trans_mat, k, m]= trans_state(isingmat, n, T, b, flip_num);
                isingmat = trans_mat;
            end
        end       
    end
    
    %将mat_list写成txt文件
    fname = sprintf('T_%.1f-num_%d-gap_%d-samlen_%d-collstep_%d.txt', T, num, gap, list_len, collect_step);
    fid = fopen(['.\',fname],'w');
    for k = 1:num
        mat = mat_list{k};
        for i = 1:n
            for j = 1:n
                fprintf(fid, '%d\t',mat(i,j));
            end
            fprintf(fid, '\n');
        end
    end
    fclose(fid);
    
    figure(1)
    scatter(steps, Clist)
    xlabel('steps');
    ylabel('C (a.u.)');
    title_string = sprintf ('specific heat curve (T = %.1f, b = %.1f)', T, b);
    title(title_string);
    
    img_step = round((num-3)/(img_num-1));
    for i = 2 : img_step : (num-1)   
        figure(i)
        ising_plot(n, mat_list{i});
        title_string = sprintf ('configuration of step %d (T = %.1f, b = %.1f)', i*gap ,T, b);
        title(title_string);
        pause(0.5);
    end
    
%     for step2stable = steps
%         for i = 1:step2stable
%             [trans_mat, k, m]= trans_state(isingmat, n, T, b, flip_num);
%             isingmat = trans_mat;
%         end
% 
%         for i = 1:list_len
%             for j = 1:collect_step
%                 [trans_mat, k, m]= trans_state(isingmat, n, T, b, flip_num);
%                 isingmat = trans_mat;
%             end
%             en = ising_energy(isingmat, n, T, b);
%             en_list(i) = en;
%         end
%     end
    
    
    
end

function E = ising_energy_in(mat, n, b)
    J = 1;
    E = - b * sum(sum(mat));
    for i = 1 : n
        for j = 1 : n
            if(i == 1)
                if(j==1)
                    sumi = mat(1, 2) + mat(2, 1) + mat(1, n) + mat(n, 1);
                elseif(j==n)
                    sumi = mat(1, n-1) + mat(2, n) + mat(n, n) + mat(1, 1);
                else
                    sumi = mat(1, j-1) + mat(1, j+1) + mat(2, j) + mat(n, j);
                end
            elseif(i == n)
                if(j==1)
                    sumi = mat(n, 2) + mat(n-1, 1) + mat(1, 1) + mat(n, n);
                elseif(j==n)
                    sumi = mat(n, n-1) + mat(n-1, n) + mat(1, n) + mat(n, 1);
                else
                    sumi = mat(n, j-1) + mat(n, j+1) + mat(n-1, j) + mat(1, j);
                end
            else
                if(j==1)
                    sumi = mat(i-1, 1) + mat(i+1, 1) + mat(i, 2) + mat(i, n);
                elseif(j==n)
                    sumi = mat(i-1, n) + mat(i+1, n) + mat(i, n-1) + mat(i, 1);
                else
                    sumi = mat(i-1, j) + mat(i+1, j) + mat(i, j-1) + mat(i, j+1);
                end
            end
            E = E - 0.5 * J * mat(i,j) * sumi;
        end
    end
    return
end